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ABSTRACT 


A description is given of the equations 
and the computer program being developed to 
model the complete Otto cycle. The program in- 
corporates such important features as: 1) heat 

transfer, 2) finite combustion rates, 3) com- 
plete chemical kinetics in the burned gas, 

4) exhaust gas recirculation, and 5) manifold 
vacuum or supercharging. Changes in thermo- 
dynamic, kinetic and transport data as well as 
model parameters can be made without reprogram- 
ming. Preliminary calculations Indicate that: 
1) chemistry and heat transfer significantly 
affect composition and performance, 2) there 
seems to be a strong Interaction among model 
parameters, and 3) a number of cycles must be 
calculated in order to obtain steady-state 
conditions . 
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THE INTERNAL COMBUSTION ENGINE has long been a 
source of frustration to those who would accu- 
rately and predictively model its behavior. It 
is not at all surprising that this should be the 
situation for the engine is inherently a very 
complex mechanical device. It is an open system 
with a movable boundary which exchanges both 
mass and energy with its environment, operates 
cyclicly over a broad range of temperatures and 
pressures, and contains matter which undergoes 
complex chemical transformations. This ubiqui- 
tous piece of hardware simultaneously embodies 
the principles of thermodynamics, chemical 
kinetics, transport phenomena, and fluid mechan- 
ics. Each separately is a challenge to the 
specialist; in concert they seem overwhelming. 
Yet, in spite of the difficulties, the modeler 
has pursued his objective impelled by a real 
need for the insight provided by model calcula- 
tions. - 

The early attempts to make quantitative 
calculations were hampered on three fronts. 

There was a serious shortage of reliable phys- 
ical data on thermodynamic, transport, and chem- 
ical rate properties, The computational capa- 
bilities were primitive by present standards. 

The numerical analysis that was necessary for 
effective computation had not yet been developed. 
To ameliorate the situation and permit at least 
rudimentary calculations, it was necessary for 
the early modelers to make gross simplifications. 
But in the Intervening years, these three ser- 
ious deficiencies have, in large measure, been 
corrected so that now 1) we; have an ample and 
growing store of physical data, 2) there exist 
readily available large-scale computational 
facilities , and 3) many of the computational 
techniques are well understood, This growth in 
capability, if anything, has been exceeded by a 
growth of interest in understanding the workings 
of the internal combustion engine. The growth 
of interest has been catalyzed largely by fed- 
eral pollution requirements and a national con- 
cern for the depletion of our petroleum re- 
sources. These factors have also been respon- 
sible for broadening the area of interest from 
only engine performance to include an interest 
in pollutant emissions . 

At this juncture it seemed appropriate for 
us to reexamine the situation with respect to 
modeling the internal combustion engine. Our 
evaluation of the problem led us to several con- 
clusions. First, we felt that any computer pro- 
gram should be independent of the particular 
physical data base used in the calculations. 

The reason for this is that although it is true 
that physical data are essential for the imple- 
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mentation of any mathematical model and do af- 
fect the results calculated by that model, they 
are not a part of the model itself but reflect 
the intrinsic properties of matter. Hence, it 
should be an easy task to incorporate additional 
or improved physical data as they become avail- 
able and, further, the program's operation 
should not be adversely affected by data revi- 
sion nor should data revision require any repro- 
gramming. All too often we have seen apparently 
insignificant data changes produce unexpectedly 
large effects on the calculated results. Our 
approach permits the user to stop worrying about 
the effect of unnecessary data assumptions built 
into the model itself. At all times, the best 
and most complete data base Is to be used in 
the calculations. The validity of this point 
of view’ is amply demonstrated by the success 
of our computer program for chemical equilib- 
rium computations (1,2)* which has achieved 
a worldwide distribution. Second, the com- 
puter program should model the complete cycle, 
that is, it should model the induction and ex- 
haust portions of the cycle as well as the 
compression and expansion phases of the cycle. 

Also it should be capable of repeatedly ex- 
ecuting the cycle so that it can converge 
to the conditions that exist during the steady- 
state operation of the cycle. Third, the com- 
puter program should be capable of executing 
a full kinetic treatment of the chemistry rather 
than using a truncated reaction mechanism, such 
as the Zeldovich mechanism, and supplementing 
this with steady-state assumptions. Steady- 
state assumptions can only be justified a_ 
pos teriorl , that is, by a comparison of the re- 
sults from a full kinetic treatment with those 
from the approximate treatment. Gelinas (3) 
made such a numerical comparison and concluded 
that the use of steady-state assumptions "is 
fraught with peril." Approximate chemistry may 
be adequate under some conditions and for some 
special purposes, as demonstrated by Annand (4) , 
but to establish this adequacy one must be able 
to do the calculation with more nearly complete 
chemistry. Of course, multicycle calculations 
and complete kinetics demand efficient computa- 
tion so as not to require an excessive amount 
of computer time. Fourth, there should be con- 
siderable flexibility in the choice of modeling Zeleznik 
functions and their parameters so that one can 
directly compare the results for different 
choices. Finally, fifth, the user should per- 3 
ceive the computer program as a black box for 
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carrying out cycle calculations Cor his choices 
oi* physical data and modeling parameters, and he 
should not have to be concerned with the com- 
plexities oC the computation. The results 
should be available not only in a convenient and 
easily digestible form, but also in a form which 
supplies as much detail as he could conceivably 
require. 

At the time we began our work, we were not, 
and still are. not, aware of any computer program 
which possessed ail of these features which we 
deem to be so desirable. Many calculations were 
based on precalculated equilibrium properties 
which were then utilized with some interpolative 
scheme instead of calculating them as needed. 

Others made wholly unnecessary assumptions about 
thermodynamic data. Still others ignored heat 
transfer effects which strongly affect the tem- 
perature and which, in turn, could have a 
drastic influence on die composition through the 
temperature dependence of the thermodynamic, 
transport, and chemical rate data. Oftentimes 
only -the closed system portion of the cycle was 
computed, based on estimated conditions at the 
start of the compression stroke. Generally, 
chemical kinetics was given short shrift either 
by ignoring it completely or else by drastically 
curtailing the number of reactions in the reac- 
tion mechanism and supplementing these with 
steady-state assumptions. Transport properties 
were totally Ignored or were represented by the 
most gross approximations. 

What are the prospects for developing a 
computer program to model the internal combus- 
tion engine which avoids those features which 
we find undesirable and incorporates those which 
we consider desirable? Clearly, we feel that it 
should be possible to accomplish the task, but 
also we realize that the foregoing considera- 
tions impose stringent constraints on the level 
of description that can be attempted in the 
modeling. Complete multicycle calculations and 
chemical kinetics separately, and certainly to- 
gether, preclude modeling with partial differ- 
ential equations because of the attendant vast 
requirements for computer time. Thus, the cal- 
culation rajst, of necessity, be based on a sys- 
tem of ordinary differential equations. This 
conclusion is reinforced by the stipulation that 

the calculations be independent of the partic- Zeleznik & McBride 

ular data base, for then not only is there an 

increase in the bookkeeping overhead associated 

with the calculation but also one cannot take ^ 

advantage of computational devices which might 

be dependent upon accidental characteristics of 

a given data base. To reap the potentially 

large benefits of broad applicability and con- 
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venience of use one must pay the price of in- 
creased computation time. This increased cost 
can be made quite modest by paying careful 
attention to numerical techniques. 

On the following pages we shall briefly 
describe our attempt to achieve the objectives 
we have been discussing. We shall also illus- 
trate the progress we have made up to this point 
with some typical calculations. However, 
neither our descriptions nor our illustrations 
should be regarded as our final words on the 
subject for we are continuing our development of 
the computer program and the analysis. Our 
ultimate goal is to use our computer program to 
predict engine behavior. To accomplish this we 
expect to use our program, together with engine 
data being generated at our laboratory, to de- 
termine suitable values for model parameters 
from experimental data as they become available. 

GOVERNING DIFFERENTIAL EQUATIONS 

For discussion purposes it is convenient to 
partition the, Otto cycle into seven segments on 
a pressure-volume indicator diagram by eight 
points, two of which coincide. This can be seen 
in Fig. 1. Each segment corresponds to a dif- 
ferent aspect of the cycle. The segment [1,2] 
represents the reduction of the residual exhaust 
gas pressure from the exhaust pressure to the 
manifold pressure. The induction of a fresh 
fuel-air-recirculated exhaust gas mixture and 
its mixing with the residual exhaust gases takes 
place on the segment [2,3], The contents of the 
cylinder are compressed by the piston along seg- 
ment [3,4] and combustion is initiated at the 
point labeled by 4. Combustion now continues on 
[4,5] where we have the coexistence of burned 
and unburned working fluids. Essentially all of 
the cylinder contents have been burned by the 
time point 5 has been reached and practically no 
unburned gas remains. The combustion products 
are then expanded along the segment [5,6]. At 
point 6 there is a sudden reduction of pressure 
to the exhaust pressure followed by the exhaust 
of the combustion products from the engine along 
segment [7,8]. A different system of differen- 
tial equations is applicable, and is integrated, 
along each segment but there are always enough 
equations to determine the thermodynamic state 
of the working fluid and its composition. The 
specification of the state demands any two inde- 
pendent thermodynamic variables. During combus- 
tion, the burned and unbumed gases each require 
two variables for a total of four. 

We shall now simply give the equations 
which are the basis for our modeling of the 
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seven segments of the cycle. The derivation and ! ; 

justification of these equations would take us j i 

too far afield and these aspects of the equa- f 

tions will not be discussed in this paper. j 

Clearly the equations should be established with 
a minimum use of restrictive assumptions to 
achieve a maximum domain of applicability, 

SEGMENT [ 1 , 2 ] AND SEGMENT [ 3 , 4 ] - Identical f 

equations are used for segment [1,2] and segment || 

[ 3 , 4 ]. The erankangle Intervals corresponding 
to these segments are 0 = 6^ <. 6 <, e 2 and 
tt » 93 £. 0 04, respectively. 


du 

1 dV Q 

(1) 

dO 

3 M [“P d0 ~ «J 

dv 

f 1 dv 1 dnl 

(2) 

d0 = 

v [v do " M dej 

dn^ 



dT~ “ 0 

SN 

T— 1 

II 

( 3 ) 


dM „ 

— — a 0 

d0 

( 4 ) 


V = V( 0 ) 

( 5 ) 


. M 2-° 

(6) 


V « ° 

( 7 ) 


SEGMENT [ 2 , 3 ] - The erankangle interval for 
this segment, which is the charge induction por- 
tion of the cycle, is 62 £• ® — 83 “ -tI * The 
relevant equations are the following. 
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V 

= V(6) 

(12) 







dM, 

de“ = 0 

( 13 ) 
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The reactants comprising the fresh fuel-air mix- 




ture are initially presumed to be at u pressure 
Pp and are iaentropicully expanded to the pres- 
sure p M prior to induction into the cylinder . 
along [2,3], 

SEGMENT [4,5] , ,'his segment is the com- 
bustion phase of the cycle and the crankangle 
interval is 64 £. 9 0^ . 


(14) 


Vlf dV l *1,, dM l’ 

i M 1 P 1 do a + ^ l 2 “ u l^ de 
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da 3 V 1 v x de ~ do 
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X » 1, 2, . . . (16) 
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dT 3 ^[“ p 2 dr"T + Ch 2 dTj^ 
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3 p 2 
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V l l + 

v 2 M 2 

= V 

(21) 


M 1 + 

m 2 - 

M 

(22) 


dM 

da 

= 0 
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M l“ 

M x (0) 


(24) 


If the products of combustion are to be in chem- 
ical equilibrium along [4,5], then this assump- Zeleznlk & McBride 
tion is used to replace the differential equa- 
tions for . 

SEGMENT [5,6] - The crank angle interval 7 
for the post-combustion phase of the cycle is 
65 S. 6 0g - 3n. The following equations are 
used to calculate the working fluid properties 
for this segment. 
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du h 1 [ dV a 

da m l P do 


(25) 

dv f l dV 1 dM 

dO “ v |_V d0 ” M dO _ 


(26) 

dn. R. 

~ X » 1, 2, . 

d0 pm ’ 

• • 

(27) 

dM _ 

sir 0 


(28) 

V = V(6) 


(29) 

M 2 e 0 


(30) 

If the products are to be in chemical equilib- 
rium alocg [5,6] then, as for [4,5], this 
assumption is used to replace the differential 
equations for n^, except for NO which Is as 
assumed to be nonreacting. 

SEGMENT [6,7] - The sudden reduction of 
the pressure to exhaust pressure is calculated 
from the following equations for 6g 3 By H 3ir. 

n X (0 7 ) " n X (t V = 0 

(31) 

s(6 7 ) - s(e 6 ) « 0 


(32) 

p( ft 7 ) - P E “ 0 


(33) 

v(6 7 )M(0 7 ) - v(0 6 )M(0 6 ) 

- 0 

(34) 

SEGMENT [7,8] - This segment 
corresponds to the exhaust stroke 
3w » 8 7 <. 0 £. flg ~ 4tt, 

of the cycle 
for which 

dh b _ 
de ffij 


(35) 

HsT 0. p(e 7 ) - 1- 

>E 

(36) 

dn.. 

d r =0 ial - 2 " 

• • 

(37) 

v(0) = V(e)/M(9) 


(38) 

v = v(e) 


(39) 


In essence, these equations imply that the 
Otto cycle is an open system only for the seg- 
ments [2,3], [6,7], and [7,8]. Chemical reac- 
tions occur only in the burned gases and there 
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only along [4,5], and [5,6]. The reaction rates 
are determined either by chemical kinetic data 
or by equilibrium thermodynamics. For the 
latter case NO is assumed to be nonreucting in 
the post combustion segment [5,6], When chem- 
ical kinetics governs the combustion process, it 
is ussutned that the unburned gas passes through 
the flame front and is injected into the burned 
gas at the equilibrium composition corresponding 
to the enthalpy and pressure of the unburned 
gas,. Heat and work can be exchanged between the 
working fluid and its environment at all points 
of the cycle. There is no heat transferred be- 
tween the burned and unburned working fluid, and 
both are at the same pressure but not the same 
temperature. 

The preceding equations are the ones we 
have selected as the governing equations for the 
seven segments of the Otto cycle, but any 
attempt to implement them in the given form 
would certainly lead to serious nuwdsricul diffi- 
culties. These equations have been converted to 
numerically useful forms} however, the details 
of the transformations would entail the discus- 
sion of a substantial amount of analysis. For 
that reason this topic will not be discussed 
here, 

MODELING FUNCTIONS - 'Hie systems of equa- 
tions we have been discussing are complete in 
the sense that their solutions can be used to 
completely characterize the Otto cycle at the 
level of description that we have chosen (ordi- 
nary differential equations). The equations 
themselves are, however, still incompletely 
specified because they contain functions , called 
modeling functions, which remain tc be defined. 

These functions include the heat loss functions 

• * • 

Q, Ql, and Qt, the mass burning function Mi, 
the volume function V, and the species rate of 
production functions R^. It is important to 
realize that while there are sometimes conven- 
tional choices which are made for some of these 
functions, it is by no means certain that the 
conventional choices are the "best" choices. It 
is certainly ooseible, and in fact likely, that 
there is no uniquely "best" choice. The ulti- 
mate criteria for a good choice are adequate 
agreement with experiment and extrapolation ade- 
quacy. Nonetheless, to perform any calculations 
one must make a choice. We must make some spe- Zeleznlk 
cific choices because we wish to present the re- 
sults of some typical calculations within the 
framework of our general approach in order to 
illustrate the versatility and computational 
costs and also to demonstrate some parametric 
effects. We do not wish to attach any particu- 
lar significance to the ad hoc choices we make 
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for these calculations , We might well make dif- 
ferent choices for subsequent calculations. The 
following equations give our choices of the 
modeling functions for this paper, 

V(e) o v c [l + (r - 1)[1 - cos 0 ; , 

+ (1 - A - t 2 sin 2 0)/e]/2| (40) 


H - 


| *» H 


e “ L/2& 

W - \) 




- a.) /* jJ} i e. state, («> 


Hjte,) + [H - MjCSfUU - - 0,»/'*,U S ( ± *s 


q - hA(T - T w ) 


Ql - hA 1 (T 1 - T w ) 
Q 2 u %<I 2 " T w ) 


A - A + 2(r - 1)V [1 - cos 

C v* 


+ (1 - AT- c 2 sin 2 6)/e]/B (46) 


“ AM^/M 

\ “ am 2 /m 


c 1 (L»/*) 1/3 <pi) 1/2 


Eichelberg (49) 


- c«(k/B) (BLtop/tip) 0 * 7 + c 3 (T 4 - T 4 )/(T - T w ) 


Annand (50) 


« c 4 (k/B)(c u/k)°’ 4 (BL UJ p/irp) 0,8 


Dittus-Boelter (51) Zeleznik & McBride 


Equation (40) for V(o) is just the usual 
volume-crankangle equation used for the four 
stroke piston engine. Equation (42) for M^(0) 
is taken from (5) . For 0 it reduces to 

the mass burning expression of Blumberg and 
Hummer (6) and for 9 ^ 0 it is the extension 




of choir formula to incorporates the relaxation 
effects of Blizard and Keck (7). The heat loss 
formulas are the conventional ones involving the 
heat-transfer coefficient. The expressions for 
Ai and A 2 were motivated by the knowledge 
that the separation of the working fluid into 
burned and unburned fractions by the flame front 
is sufficiently indistinct that and A 2 

canno* be determined precisely. Further, they 
possess the desirable property that as the mass 
approaches zero, so does the Sieat transfer sur- 
face, We assume the usual rate expressions of 
chemical kinetics for %. If C]_ » 1. 25x10-5, 
c 2 in the interval (0.35, 0.8), 
c-. « 1.025«10“ 13 , and c, - 0.023, then the for- 
mulas for the heat transfer coefficient, in 
cal-citr^-sec“^“K.“^, are the correlations of 
Elchelberg (8), Annand (9), und Dittus-Boelter. 

The latter is a correlation for heat transfer 
accompanied by the turbulent flow of a fluid in 
a pipe. Each of the heat transfer correlations 
has features not possessed by the other two. 

The Elchelberg form is the one most commonly 
used and has no composition dependence. The 
Annand form has a radiation term and depends on 
composition through the thermal conductivity. 

The Dittus-Boelter form depends on composition 
through both the thermal conductivity and the 
shear viscosity. 

To conclude the discussion of modeling 
functions, we wish to point out that many model- 
ing functions, like the formulas for the heat 
transfer coefficient, often contain parameters 
which are determined from experimental data by a 
comparison with calculated results. But since 
the calculated results can be strongly affected 
by the governing equations and by the physical 
data, one should realize that the parameters may 
have to be redetermined as the ability to calcu- 
late improves. That is, experimentally deter- 
mined parameters may be model sensitive. 

PHYSICAL DATA 

An extremely important part of any modeling 
computer program is the physical data base which 
is used in the calculations, To model the in- 
ternal combustion engine with tiie equations of 
the previous section we require thermodynamic 
data for each of the chemical species which Zeleznik 

appear in the working fluid. For each species, 
we must know the heat capacity, the enthalpy, 
and the entropy. All of the necessary thermo- 11 
dynamic computations are performed with our own 
equilibrium computations program (2). Hence, 
the thermodynamic data are just those which we 
supply with that program and these are supple- 
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men ted by some additional duta for the fuels 
which ate needed for the calculations’ to be 
reported In this paper. 

The transport properties of the working 
fluid appear in two of the three formulas for 
heat transfer coefficients. The shear viscosity 
and thermal conductivity data for the pure spe- 
cies were based on experimental measurements 
where possible and, when necessary, this vjuh 
augmented with theoretically estimated proper- 
ties, Except for the thermal conductivity of 
carbon dioxide, the transport properties came 
from (10 to 13).. The thermal conductivity of 
carbon dioxide was obtained from a smoothed 
curve through experimental data from 16 refer- 
ences published from 1950 to 1968. The mixture 
properties were calculated from the data for the 
pure species with the mixture formula for vis- 
cosity suggested by Bromley and Wilke (14) and 
the mixture formula of Lindsay and Bromley (15) 
for thermal conductivity. The mixture thermal 
conductivity formula was modestly altered so us 
to make use of information contained in the 
shear viscosity data. 

The chemical kinetics data consisted of the 
rate constants and equilibrium constants for the 
19 forward reactions listed in Table 1. This 
table contains only homogeneous reactions and we 
will neglect-, all mil inactions » It should be 
noted that four of the listed reactions involve 
third bodies and so each of these reactions cor- 
responds to 14 forward reactions, one for each 
possible third body. The rate constant datu for 
the individual forward reactions came from (16) 
for all reactions except reaction 13 and reac- 
tion 19. The rate constants for these two reac- 
tions came from (17) and (18), respectively. 

The equilibrium constant data used to compute 
the rate constants for the backward reactions 
from the forward rates came from the same ther- 
modynamic data mentioned previously. The net 
production rates were calculated by the 

method described by Blttkcr and Scullin (19) and 
we used the appropriate portions of their com- 
puter program to carry out the calculations. 

SOME CALCULATED RESULTS AND THEIR IMPLICATIONS 

THE STANDARD PROBLEM - Our purpose in pre- 
senting the results of some numerical calcula- Zeleznik & McBride 
tiona is 1) to display the capabilities of the 
present version of our Otto cycle computer pro- 
gram, 2) to make an assessment of the validity 
of equilibrium thermodynamics during combustion 
as an approximation for chemical reaction 
effects, 3) to try to make a judgment about the 
degree of interaction of model parameters, 4) to 


determine Che sequence and number of cycles 
which should be calculated to determine steady 
state conditions and 5) to make some estimate 
of the computation times involved. We do not 
intend to carry out a detailed parametric study 
at this time because we feel that too many of 
the model parameters given in the literature are 
so poorly estimated that such a calculation 
would be premature. 

Our computer program for modeling the Otto 
cycle is already very versatile in its present 
stage of development. It can accommodate so 
many parameters that it would be a herculean 
task to illustrate the effects of each one. 

Since our purpose is only to make some general 
observations about the calculation^ , we shall 
arbitrarily fix a number of parameters at what 
might be considered reasonable values. The 
assigned values are shovm in Table 2. These 
values, together with Pp “ Pe “ nun 
Pjj - 250 mm Hg, a stoichiometric air-fuel ratio, 
no recirculated exhaust gas, t ° 0, and a fre- 
quency of 500 revolutions per minute will define 
a standard set of conditions for all comparison 
purposes. The heat transfer parameters eg and 
c, were chosen so that the heat loss was iden- 
tical for all three heat transfer correlations 
for the standard conditions and a particular 
cycle in the sequence of cycle calculations. 

The discussion of the calculated results 
will be facilitated if we first introduce some 
terminology associated with the sequence of 
cycle calculations executed by the computer pro- 
gram. To reach steady state conditions with the 
least amount of computer time usage, it is de- 
sirable to begin the calculation with the most 
highly idealised model and sequentially intro- 
duce departures from the ideal. The greater the 
idealization the simpler the calculation and, 
hence, the shorter the computation time. Thus, 
the computation begins with a set of n-j- ideal 
cycle calculations corresponding to equilibrium 
combustion at constant volume and zero heat 
transfer during the cycle. The actual number jf 
ideal cycles computed is the number required to 
achieve steady state operation, defined by the 
condition that thii cycle to cycle change in 
exhaust gas temperature at point 7, T(O^) , is 
less than 1 K. These ideal cycles are then fol- 
lowed by ng cycles which are somewhat less Zeleznik 

Idealized and include the effects of heat trans- 
fer and a finite combustion rate. The composi- 
tion of the burned gas is still governed by 13 

equilibrium thermodynamics. These cycles in 
turn are followed by ng cycles where chemical 
kinetics is used to determine the burned gas 
composition in place of equilibrium thermo- 
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dynamics . The nonnegative integers n E and n^ 
are arbitrary and are limited only by the avail- 
able computation time. A given calculation can 
then be labeled by the triplet of integers (n^, 
n E> • Th® heat transfer parameters C£ and 
04 were chosen so as to give identical heat 
losses for the calculation (nj, 3, 0) of the 
standard problem, that is, at the conclusion of 
the third equilibrium cycle. 

Some calculated results for the standard 
problem, using the Eichelberg form of the heat 
transfer correlation, are shown xn graphical 
form as Figs. 2 to 12, These figures are com- 
puter generated plots which are pretty much 
self-explanatory and are the results of the 
(3, 3, 0) calculations. That is, these results 
are for the third equilibrium combustion cycle 
which was preceded by three ideal cycles. Each 
figure is labeled by the date of the calcula- 
tion, a case number and a cycle number. The 
first digit of the case number identifies the 
form Of the heat transfer correlation (1 « 
Eichelberg, 2 « Dlttus-Boelter, 3 ** Annand) and 
the next two digits give the revolutions per 
minute in hundreds of revolutions. The cycle 
number Is the number of cycles calculated after 
the steady state ideal cycle. The term "cyl- 
inder average" used in Figs. 7 and 8 refers to 
an average over the burned and unburued gases. 

From these figures, one can see that the 
values of all properties agree at 0=0° and 
0 ** 720° and hence, these figures correspond to 
a steady state where combustion is controlled by 
thermodynamics. The curves in Figs. 2, 3, 4, 
and 9 possess discontinuities at 0 - 540° 
which occur because of the sudden reduction of 
the pressure to pg. Figures 9 to 11 each have 
two curves. One curve is either the heat or the 
work and the other is 100 times its rate ct 
change with respect to crankangle. On the 
boundaries of the Otto cycle segments, that is, 
at the distinguished points of Fig. 1, the rate 
of change is the actual derivative at the ini- 
tial point of the segment. At all other points 
the rate of change is the mean rate of change 
obtained by averaging the rate of change over 
intervals determined by the computer program. 

The discontinuity in the heat loss rate at the 
end of combustion (point 5) arises for this 
reason and occurs because of the rapid change 
in the heat loss rate in this portion of the 
cycle. 

CHEMICAL KINETICS EFFECTS - What is the 
effect of replacing equilibrium thermodynamics 
with chemical kinetics and calculating two addi- 
tional cycles for the same problem, that is, do- 
ing the calculation (3, 3, 2)? For this par- 


ticular calculation, using the Eichelberg heat 
transfer correlation, there 1b essentially no 
change in temperature, pressure, work, and heat 
loss. There are, however, changes in the con- 
centrations of the pollutants carbon monoxide 
and nitric oxide. These are displayed In 
Figs. 13 to 16 for the calculation (3, 3, 2), 
Again, because the values at 0 » 0° and 
0 = 72 0° are equal, these results correspond to 
a new steady-state where combustion is con- 
trolled by chemical kinetics, A comparison of 
these figures with those given as FigB . 5 to 8 
snows that a considerable difference exists be- 
tween pollutant concentrations calculated on the 
basis of thermodynamic equilibrium and Chose 
calculated on the basis of chemical kinetics. 
Chemical kinetics calculates higher concentra- 
tions for both pollutants. The differences in 
results between these two assumptions are fur- 
ther emphasized by the data given in Table 3, 

It shows the dependence of the calculated re- 
sults on the cycle number and the heat transfer 
correlation for the standard problem. Here we 
again see evidence that a steady state has been 
achieved in the standard problem for (3, 3, 0) 
for the four kinds of heat loss calculations and 
all three heat transfer correlations do indeed 
give the same heat loss on the third equilibrium 
combustion cycle. Two additional kinetic com- 
bustion cycles suffice for the attainment of a 
new steady-state for the quantities listed in 
Table 3. It is obvious from the data in Table 3 
that heat losses cannot be ignored because they 
profoundly alter the calculated results as can 
be seen by comparing the zero heat transfer cal- 
culation with the other three calculations. 

A new and somewhat surprising feature 
appears in Table 3. While there is very little 
interaction between the chemistry and the 
Eichelberg heat transfer correlation, chemistry 
does have an appreciable effect on the heat 
losses calculated by the other two heat trans- 
fer correlations. Also changes were produced In 
the indicated mean effective pressure, the indi- 
cated specific fuel consumption as well as the 
temperatures and pressures. This points up the 
model sensitivity of heat transfer parameters 
and should cause one to question the correctness 
of heat transfer parameters determined from the 
combination of experimental data and model cal- 
culations which ignore chemical kinetics. But, 
to the best of our knowledge, all calculations 
prior to those reported here have not had the 
combination of a general heat transfer correla- 
tion with complete chemical kinetics. Thus it 
follows that those correlations that already 
exist in the literature must be used with con- 
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siderable caution. Further, the strength of. the 
interaction between chemistry and heat transfer 
should cause one to suspect, and be wary of, 
strong interactions among other parameters. For 
this reason , it seems probable that very few 
model parameters can be isolated for experi- 
mental determination. It is more likely that 
most: model parameters must be evaluated simul- 
taneously • from all of the available experimental 
data. 

Table 3 leads us to two other Important 
conclusions. First, the (n^, 3, 2) calculation 
seems to be adequate for obtaining an estimate 
of the true steady-state conditions. Second, 
the estimation of the concentration of the pol- 
lutants CO and N0 X by equilibrium combustion and 
frozen NO in the post combustion region may be 
in serious error. This latter point is corrob- 
orated by the data in Table 4 which shows the 
amount of N0 X , expressed as grams of N0 X per 
kilogram of burned gas, calculated by the two 
methods discussed in this paper as well as the 
values calculated if the NO were frozen at the 
peak temperature. The first part of the table 
shows the dependence on frequency and the choice 
of heat transfer correlation while the second 
part shows the dependence on equivalence ratio 
for the Eichelberg heat transfer correlation at 
a fixed frequency of 3500 revolutions per min- 
ute, From Table 4 we see that freezing NO at 
the end of combustion is superior to freezing it 
at the peak temperature as a device for estimat- 
ing N0 X ; however, neither technique gives a par- 
ticularly good approximation to the kinetically 
calculated values. 

Exhaust gas recirculation has been used as 
a method for reducing the emission of N0 X but it 
should, and does, affect the concentrations of 
other species in the exhaust. For case 135 the 
specific productions of NOx and CO in g/kw-sec 
are 4.62*10" 3 and 3.13xl0“3 } respectively, for 
(3, 3, 0) and 4.40xl0" 3 and 3.12x10"-* for 
(3, 3, 2). The addition of a sufficient amount 
of recirculated exhaust gas so that it repre- 
sents 15% of the working fluid by mass (case 
13515) changes these values to 2.53x10-3 and 
1.57 *10 *“ 3 for (3, 3, 0) and 1.93x10"-* and 
1.68X10" 3 for (3, 3, 2). Thus for both equi- 
librium combustion and kinetic combustion the 
recirculation of exhaust gas significantly re- 
duces the production of N0 X and CO. 

One of the most difficult pollutants to 
handle in a model of the Internal combustion 
engine is the unbumed hydrocarbons, Zeleznik 
(5) suggested that the relaxation approach of 
Blizard and Keck (7) might offer an approach to 
the problem. To see if unburned hydrocarbons 


Zeleznik & McBride 


16 


iukjJi 


could be made to appear In the exhuust we re- 
peated chfie 135 but set t » 0.0009. The result 
was that for this case, designated 13509, 

1.9 parts per million hydrocarbon appeared in 
the exhaust. Hence, this is indeed a means for 
introducing hydrocarbon into the exhaust. The 
only question is whether or not one can deter- 
mine the parameter x ao as to reproduce exper- 
imental unbumed hydrocarbon data. This can 
only be answered after an extensive analysis of 
calculated results and experimental data. 

The concentrations of CO and NO generally 
behave qualitatively as shown in Figs , 5 to 8 
and Figs, 13 to 16; there are, of course, sig- 
nificant quantitative differences. But we have 
also seen cases where the qualitative behavior 
differs appreciably from what we have shown. As 
an example of a different kind of behavior, we 
show in Figs. 17 and 18 the concentrations of CO 
and NO in the burned gas of the calculation 
(3, 3, 2) of case 13509 and case 13515, respec- 
tively. Here the carbon monoxide decreases ini- 
tially before it begins its rise to a maximum. 

The nitric oxide rises to, and then remains at, 
its maximum value, 

DETERMINATION OF HEAT TRANSFER PARAMETERS - 
What are the prospects for determining param- 
eters in heat transfer correlations from experi- 
mental data? The prospects would be excellent 
if the calculated results were strongly depend- 
ent upon these model parameters through some 
externally controllable parameters such as fre- 
quency and air-fuel ratio. We have already seen 
from Tables 3 and 4 that the composition is mod- 
erately dependent on both the form of the heat 
transfer correlation and the frequency and is 
strongly dependent on equivalence ratio. Hence 
data on pollutant concentrations in the exhaust 
might be of some help in evaluating heat trans- 
fer parameters . What about the frequency de- 
pendence of heat loss and indicated work and its 
relationship to the heat transfer correlation? 

This is shown in the first part of Table 5 where 
we give the frequency dependence of the heat 
loss and the indicated work for the three corre- 
lations used in this paper at the steady-states 
for equilibrium and kinetic combustion. From 
these data we see that the frequency dependence 
of the heat loss and work is greatest when the 
calculation Is made with the Eichelberg heat Zeleznik 

transfer correlation and weakest for the calcu- 
lations made with the Dittus-Boelter correla- 
tion. But in any event, the frequency depend- jj 
ence of the calculated results is sufficiently 
different for the three correlations that it 
should be possible to use experimental heat loss 
and work data to choose from among various cor- 




McBride 


relations. The information derived front the 
frequency dependence could be augmented by ex- 
ploring the dependence of the heat loss and work 
on the air-fuel ratio. The kind of variation 
predicted with the Eichelberg correlation is 
shown in the second part of Table 5 and it seems 
reasonable that the results calculated with the 
other correlations should be affected even more 
strongly because of their greater sensitivity 
to the chemistry. 

The Annand form for the heat transfer coef- 
ficient is the only correlation which tries to 
make an allowance for radiation. Is this T^ 
term likely to contribute significantly to the 
heat loss? Since heat loss per cycle decreases 
as the revolutions per minute increases, the gas 
temperatures should be the highest at 5000 rpm. 

Hence, the radiation contribution should 
approach its maximum value under these condi- 
tions. We, therefore, repeated the calculation 
for case 350 but eliminated radiation effects by 
setting the constant c^ to zero. This de- 
creased the heat loss by about 1.5% and in- 
creased the work by about 0.3% for both kinetic 
and equilibrium combustion chemistry. Thus, the 
effect of radiation Seems to be quite small and 
probably cannot be determined from experimental 
data with any real accuracy. 

COMPUTATION TIMES - It is very difficult to 
make any precise statements about computation 
times since they are so variable not only from 
one machine to another but also from one instal- 
lation of the same machine to another. All our 
computations were carried out on an IBM 350/67 
duplex with the TSS/360 operating system. This 
computer has virtually unlimited storage 
capacity but is generally considered to be a 
relatively slow machine by current standards. 

The computation times we will cite are for an 
unoptimized program used on this machine and we 
could have reduced the times by about 10 to 
15% with relatively little effort. 

The ideal cycles typically take about 
2.5 seconds of computer time. Cycles in which 
combustion is computed by equilibrium thermo- 
dynamics take about 1 to 2 minutes of computer 
time. Those computed with the Eichelberg heat 
transfer correlation are closer to one minute 
while the other two correlations require times 
closer to 2 minutes. The cycles for which the Zeleznik 
composition is determined by chemical kinetics 
exhibit the greatest variability in the usage of 
computer time and also somewhat longer computa- 18 
tlon times. These times span the range from 
about 1.5 minutes to about 4 minutes. For com- 
parable conditions the calculations using the 
Eichelberg heat transfer correlation require the 
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lease time. Generally, the computation time 
decreases as the frequency increases and it 
tends to decrease as one moves away from stoi- 
chiometric conditions. 

It is gratifying t0 find that the computa- 
tion times for a full kinetic cycle calculation 
are not drastically different from the computa- 
tion times for an equilibrium cycle calculation. 

The increased computation cost is quite modest. 

In fact, the cost is so modest that it becomes 
difficult to justify the use of approximate 
kinetics, such as the Zeldovich mechanism 
coupled with steady-state assumptions, because 
the sole purpose of these approximations is to 
simplify the computation and thus reduce the 
computation times. The increased coat is cer- 
tainly warranted since it dispels any doubts >/ 

one would have about results based on the use ji 

of steady-state assumptions. /j 

CONCLUSIONS if 

Based on our calculations it is possible 
to draw a number of conclusions about modeling 
the complete Internal combustion engine cycle. 

First, and most important, it is feasible and 
practical to model the complete cycle with 
ordinary differential equations with both ex- 
panded finite rate chemistry during the com- 
bustion and post-combustion phases and heat 
transfer. Second, chemistry and heat transfer 
significantly affect not only composi.tipn but 
also measures of cycle performance s,uch as 
work. Third, there appears to be a sufficiently 
strong Interaction among model parameters so as 
to preclude the isolation of Individual param- 
eters for experimental determination. It may 
be necessary to evaluate a number of parameters 
simultaneously. Fourth, in the sequence of 
cycle calculations discussed here, it seems that 
one can adequately determine steady-state con- 
ditions with the calculation (nj,3,2). Fifth, 
the estimation of pollutant concentrations by 
neglecting finite rate chemistry or heat trans- 
fer is questionable. 

ACKNOWLEDGEMENT 

We wish to thank Helen Hu for writing the 
computer plotting routines that we used to gen- Eeleznik £ McBride 
erate Figs. 2 to 18. 


19 


SYMBOLS 


A 


B 


li 

h ( 

IT 

k 

L 

£ 

M 


U A 

n x (0) 



P 

PE 

P F 

Pm 

Q 


surface area for heat transfer 
surface area of combustion chamber 
bore 

heat capacity at constant pressure 

enthalpy per unit mass 

enthalpy per unit mass for fresh charge 

heat transfer coefficient (gas side) 

thermal conductivity 

stroke 

connecting rod length 

mass of the working fluid 

moles of species \ per unit mass 

moles species A per unit mass for 
fresh charge 

equilibrium composition from combustion 
of unburned gas at ^2 *^2 

pressure 

exhaust pressure 
fuel pressure 
manifold pressure 
heat loss rate 


R X 

r 


s 


T 


T 


w 


u 


volumetric production rate of species 
compression ratio 
entropy per unit mass 

temperature Zeleznik 

effective wall temperature for heat 
trans for 

20 

internal energy pet unit mass 


a McBride 


V 


volume 


volume of combustion chamber 
v volume pet unit mass 

y crankangle 

i “ 1, 2, ... 8 crank angles corre- 
sponding to the 8 distinguished points 
on the indicator diagram given in Fig. 1 

0 t mi 

Oj crankangle for which all unburned gas 

has passed through the flame front 

y shear viscosity 

p mass density 

t relaxation time for combustion 

w angular frequency 

SUBSCRIPTS AND SUPERSCRIPTS 


1 burned gas 

2 unburned gas 
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Table 1 ~ Summary of Reaction Mechanism and Rate Constants 
k for Klnetically Controlled Combustion 

„ , Third k « AT n e cm 3 /mole~soc 

Reaction body, 

number Reaction M A n E 


1 

N + OH -y NO + H 

4*1C 13 

0 

0 

! 2 

M + 0 + 0 > 0 2 + M Q, 

1.38*10 18 

-1 

340 

3 

N + NO + N 2 + 0 

1.6*10 13 

0 

0 

; 4 

N 4* 0 2 *+■ NO + 0 

6.6*10 y 

1 

6 260 

s 

0 + NO, NO + 0 2 

1*10 13 

0 

596 

fa 

' 

M + NO + 0 **■ N0 2 + M Q, 

1.1*10 15 

0 

-1 870 

7 

H + 0 2 + OH + 0 

zmqU 

0 

16 800 

' > 8 

0 + U 2 -*• OH + 11 

1,8*10 10 

1 

8 900 

9 

H + Ii0 2 ■> OH + Oil 

2.5X10 14 

0 

1 890 

v ID 

OH + H0 2 * H 2 0 + 0, 

5«10l3 

0 

994 

11 

o + ho 2 "- > - OH + o 2 

4.8xio13 

0 

994 

j: 12 

M + H + 0 2 ** H0 2 *"+ M Ar 

2.4*10 15 

0 

-5 76 

13 

H0 2 + NO •+ OH + NO, 

1.2*10 13 

0 

2 390 

! 14 

H + N0 2 ■* NO + OH 

3.5x101-4 

0 

1 470 

i !■ 15 

0 + ll 2 0 OH + OH 

6.6*10 13 

0 

18 370 

; 16 

Hy + OH * H,0 + 11 

2.2«10l3 

0 

5 150 

ii 1? 

CO 4- OH ■* C0 2 + H 

3.1X10 11 

0 

596 

i : 18 

M+H+ OK *>■ H 2 0 + M Ar 

1.4*10 22 

-2 

0 

I; 19 

o + cq 2 ■> CO + o 2 

1.9xiq13 

0 

54 150 

r 

Third body efficiencies are 1.0 

except for the 

following: 

j 

Reaction M = 0 2 M = l^O M 

“ C0 2 H = N 2 

M 

= CO 

]i 

12 2 32.5 

7.5 2 

2 


j! 

18 1.6 20.0 

■'.5 1.6 

1 

.6 
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Table 2 - Some Parameter Alignments 


Fuel composition, 
mass percent 


n-Octune 42 
Toluene 38 
Benzene 18 
1-Octene 2 


Combustion 

parameters 

Q4 « 350° 

0 £ - 0/ « 30° 
Ml(e 5 )/M(ir) « 0.9999 
T « 0 


Air composition, Heat transfer 


mole percent 


parameters 


Nitrogen 

78.0881 

C 1 

» 1.2504xl0~ 5 

j| Oxygen 

20.9495 

c 2 

** 0.3872 . 

« 1.0253xlO" XJ '! 

j'j Argon 

0.9324 

c 3 

]| Carbon 

j! dioxide 

0.03 

c 4 

« 0.73 f, [, 

fl 

i . 


Engine and 
operating parameters 

B - 10.922 cm (4.3 in.) 

I “ 10.312 cm (4.06 in.) 

A => 17.145 cm (6.75 in.) 
r “ 10-5 

A c = 129.032 cm 2 (20.0 in. 2 ) 

p(02> “ Pm 
T w - 360 K 
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Table 3 - The Effect of Heat Transfer on Calculated Results 


1 * 
le 

Ave rage 
exhaust 
temp. , 

K 

KO x 

at 500 

CO, 

pprt 

Rev/min 

Hear loss, 
Joules 

IMEP , 
atm 

ISFC, 

g/ktt-BCC 



Hero Heat 

Transfer 




1 

1229.4 

11 640 

0 

0 

6 

.4898 

0.04867 

2 

1470.7 

8 210 

2380 

0 

5 

.9411 

.05152 

3 

1469.9 

7 880 

6170 

0 

5 

.7902 

.05210 

4 

1472.4 

6 159 

5012 

0 

5 

.8042 

.05135 

3 

1471.9 

6 150 

5004 

0 

5 

.8062 

.05137 

6 

1471.8 

6 ISO 

5004 

0 

5 

.8063 

.05137 

7 

1483. S 

3 833 

3266 

0 

5 

.8320 

.05114 

8 

1485.7 

3 835 

3264 

0 

5 

,8260 

.05117 


Eichelberg Heat Transfer Correlation 


4 

798.0 

4 

197 

1800 

604.83 

4.8163 

0.06521 

5 

783.6 

3 

199 

1288 

560.11 

4,7310 

.06392 

6 

783.1 

3 

170 

1274 

558.89 

4.7289 

.06390 

7 

781.6 

3 

460 

2376 

558.38 

4,7182 

.06405 

8 

781.9 

3 

461 

2372 

558.24 

4.7197 

.06403 


Dittus* 

-Boelter 

Heat 

Transfer Correlation 


4 

798.0 

4 

053 

1788 

593.37 

4.5843 

0.06600 

5 

772.9 

3 

000 

1202 

560.68 

4.5479 

.06478 

6 

771.6 

2 

950 

1174 

558,88 

4.5426 

.06475 

7 

724.8 

3 

453 

2380 

590.24 

4.A621 

.06592 

8 

723,1 

3 

401 

2347 

588.12 

4.4576 

.06583 


Annand Heat Transfer Correlation 


4 

784.5 

4 

340 

1870 

599.75 

4.7386 

0.06487 

5 

763.6 

3 

205 

1264 

560.60 

4.6685 

.06365 

6 

762.7 

3 

160 

1241 

558.94 

4.6637 

.06364 

7 

683.2 

3 

476 

2416 

625.25 

4.4300 

.06699 

8 

681.3 

3 

388 

2352 

620.72 

4.4235 

.06682 


Cycles 1-3 Instantaneous burning and equilibrium chemistry 
Cycles 4-6 Finite rate burning and equilibrium chemistry 
Cycles 7-8 Finite rate burning and finite rate Zeleznlk 

chemistry 
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Rev/min 

500 

1300 

3500 

3000 


Table 4 - Predicted NO x (g NO x /kg burned gas) by Different Methods 


Kinetic combustion 

Equilibrium combustion 

Peak 

tempo ruture 

Eiehel- 

Dlttus- 

Eiehel- 

Dittus- 


Elciiel- 

Dittus- 


berg 

Boulter Annand 

berg 

Boelter 

Annand 

berg 

Boelter 

Annand 

3.55 

3.49 3.48 

3.25 

3.03 

3.24 

4 .65 

4.53 

4.60 

4.45 

3.94 4.08 

4.51 

3.54 

4.02 

5. 55 

4.86 

5.21 

4.91 

4.11 4.35 

5.68 

3,89 

4.49 

6.67 

5.20 

5.53 

5.02 

4.10 4.43 

5.3? 

4.02 

4.65 

6.24 

5.31 

5.71 

Eichelberg Heat Transfer 

Correlation, 3500 

Revolutions per 

Minute 



Equivalence Kinetic Equilibrium Peak 

ratio combustion combustion temperature 


0.6 

7.50 

.9 

7.11 

1.0 

4,91 

1.1 

2.29 

1.2 

.86 


7.97 

9.12 

7.43 

8.41 

5.68 

6.67 

2.45 

3.12 

1.17 

1.79 
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Table 5 - The Interaction of Chemistry and lleat Transfer 


A. frequency Dependence of Heat Loss and Indicated Work in Joules 


Equilibrium combustion Kinetic combustion 


Rev /min 

Eichelberg 

DittuS“Boelter 

Annand Eichelberg 

Dittus-Boelter Annand 




Heat loss 




500 

558.89 

558.88 


558.94 

558.24 

588.12 

620.72 

1500 

366.15 

491.12 


447.72 

366.13 

520.15 

512.93 

3500 

243.27 

439.60 


369.83 

243.15 

467.25 

431.89 

5000 

201.37 

418.44 


339. Sd 

20Q.66 

444.73 

399.62 




Indicated work; 

'.if 

A, 


500 

462.92 

444.69 


456.54 

462.02 

436.36 

433.03 

1500 

524.50 

472.39 


497.95 

524.02 

465.74 

478.90 

3500 

548.68 

489.97 


519.70 

547.80 

484.30 

503.86 

5000 

554.77 

496.52 


526.74 

553.64 

491.12 

512.02 

B. Equivalence Ratio Dependence of Heat Loss 

i and Indicated Work in . 

Joules 


for the Eichelberg Heat 

Transfer Correlation at 3500 Rev/Mln 



Equivalence 

Equilib rium 


Kinetic 



ratio 

combustion 


combus tlon 





Heat Loss 





0.8 


185.04 


185.09 



• 

9 


216.66 


216.66 



1 . 

0 


243.27 


243.15 



1 . 

1 


233 . 01 


233.11 



1 . 

2 


220.67 


219.95 





Indicated work 





0 . 

8 


468.11 


468.54 




9 


512.49 


512.60 



1 . 

0 


548.68 


547.80 



1 . 

1 


550.73 


551.33 



1 . 

2 


547.17 


547.65 
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Fig. 3 - Calculated cylinder pressure 
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Fig, 8 - Calculated nitric oxide concentration 
in the cylinder for equilibrium chemistry 
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Fig. 11 - Calculated pump work and pump work 
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Fig. 14 - Calculated burned gas nitric oxide 
concentration for finite rate chemistry 
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Fig. 15 - Calculated carbon monoxide concentra- 
tion in the cylinder for finite rate chemistry 
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